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Modifying  Sensitivity/Specificity  for  Sensors  Using 
Positive  and  Negative  Predictive  Power  Measures 
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Abstract  —  In  the  collection  of  data  from  sensors  in  the  field, 
the  u  ncertainty  in  th  e  d  ata  may  comp  romise  th  e  ab  ility  to 
accurately  p  redict  th  e  s  tate  of  a  s  ystem.  H  erein  th  e  s  tandard 
signal  detection  theory  problem  is  examined  when  nonstationary 
effects  may  occur  in  the  data  from  the  sensors.  The  use  of  PPP 
(positive  predictive  power)  and  NPP  (negative  predictive  power) 
adds  a  ne  w  view  point  on  how  to  modify  se  nsitivity  an  d 
specificity  mea  sures  in  d  ecision  making  involving  multiple 
sensors.  Th  is  is  es  pecially  t  rue  when  s  tationary  p  roperties  in 
received  data  may  be  violated. 

I.  Introduction 

In  standard  signal  detection  theory  (SDT)  analysis,  the 
presumption  is  usually  that  the  process  acts  in  a  stationary 
manner.  The  noise  processes  which  permeate  the  data  may  be 
presumed  to  be  Gaussian  with  a  constant  mean  and  variance. 

In  real  applications,  however,  the  actual  uncertainty  modeled 
via  randomness  may  act  more  like  a  stochastic  process  with  a 
time  varying  mean  and  standard  deviation.  These  variations  in 
the  first  and  second  statistical  moments  commonly  occur,  e.g. 
when  sensors  heat  up,  cool  down,  or  have  other  time 
variations.  The  sensors  may  even  fail.  If  the  quality  of  data 
are  changing  from  the  sensors,  one  must  modify  the  decision 
making  process  to  account  for  uncertainty  induced  by  a 
compromised  time  varying  data  stream.  This  paper  will 
examine  decision  making  when  the  quality  of  data  are 
adapting  with  time  and  use  measures  such  as  PPP  or  NPP  to 
improve  upon  current  decision-making  methodologies  based 
on  static  assumptions. 

II.  Background  material 

(2. A)  Definition  of  a  Stochastic  Process 

Figure  (1)  portrays  the  concept  of  a  stochastic  process  [1]. 

The  x  axis  is  the  independent  variable,  time.  The  vertical  axis 
is  for  the  dependent  variable  of  the  stochastic  process 


M.  A.  Vidulich,  V.  S.  Finomore 
711  HPW  AFRL  WPAFB 
Dayton,  Ohio  USA 

y(t).  This  stochastic  process  may  have  a  first  moment  (mean 
p(t))  which  changes  with  time.  The  square  root  of  the  second 

moment  (standard  deviation  aft))  also  may  vary  with  time. 
These  variations  are  portrayed  in  Figure  (1)  on  the  vertical 
axis  in  the  form  of  probability  density  functions  which 
changes  temporally. 


■A.  Stochastic  Process 
Figure  (1)  -  A  Rendering  of  a  Stochastic  Process 
The  basics  of  signal  detection  theory  will  be  discussed  next. 
(2.B)  -  A  Signal  Detection  Theory’  Framework 


Ground  Truth  =  True  State  of  the  World 


Yes  ' 

No 

Test 

Yes 

H  =  hits 

FA  = 

False  Alarms 

Result- 

No 

M  = 

Misses 

CR  = 

Correct 

Rejects 

Figure  (2)  —  Confusion  Matrix  for  Statistical  Test 
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Figure  (2)  portrays  a  confusion  matrix  for  a  binary  detection 
task.  In  Figure  (2),  let  H  denote  the  number  of  hits,  CR  is 
correct  rejects,  FA  denotes  false  alarms  and  M  represents  the 
number  of  misses.  The  standard  definitions  that  will  be  used 
here  are  described  in  equations  (1-  4)  [2]: 


Sensitivity  Sn 

H  +  M 

(1) 

CR 

Specificity  ~Sp~ 

CR  + FA 

(2) 

H 

Positive  Predictive  Power  =  PPP  =  - 

Ha  FA 

(3) 

CR 

Negative  Predictive  Power  =  NPP  = - 

CR  +  M 

(4) 

(Note  Sn  ±  PPP  and  Sp  A  NPP  from  equations  (1-4).) 

(2.  C)  The  Use  of  PPP  and  NPP  to  Assist  in  Decision  Making 
Figure  (3)  shows  the  typical  statistical  test  for  a  binary 
detection  task  and  the  appropriate  parameters  to  discern 


The  Likelihood  Ratios  =  h 

Figure  (3)  —  Classical  Binary  Decision  Making  with  Optimality 

objects  with  the  signal  to  noise  ratio  d’  and  bias  parameter  c 
as  noted.  In  the  medical  diagnostics  area  physicians  now 
highly  embrace  the  concepts  of  PPP  and  NPP  rather  than  just 
using  d’  and  c  for  several  important  reasons:  (1)  The  data 
properties  may  change  with  time  with  d’  and  c  adapting  [3,4]. 
Also,  (2)  d’  and  c  are  variables  calculated  with  a  static 
analysis,  however,  the  underlying  system  may  experience 
dynamic  changes  requiring  modification  of  the  decision 
making  process.  Evidence  of  nonstationary  behavior  in 
decision  making  is  pervasive  in  the  literature  on  experimental 
psychology. 


time,  but  PPP  and  NPP  may  remain  constant.  Also  in  [6,  7] 
evidence  suggests  that  dynamic  changes  occur  in  decision 
making,  especially  in  studies  in  vigilance  and  fatigue.  In  [4]  it 
was  shown  that  PPP  may  be  constant  (performance  feedback 
provided)  but  NPP  may  decrease  with  time.  Thus  PPP  more 
specifically  delineates  what  types  of  decisions  have  higher 
credence.  This  is  a  stronger  result  than  what  is  provided  by 
the  typical  static  measures  from  signal  detection  theory  and 
has  achieved  high  acceptance  in  the  medical  community  [10]. 

III.  Some  Key  classical  relationships 

To  relate  the  terms  in  equations  (1-4)  to  some  of  the  well 
known  quantities  in  statistical  testing  for  binary  hypotheses, 
the  elements  of  figure  (2)  can  be  described  via  [8]: 

Fj4.  /r\ 

a  = - =  1-  specificity  P) 

CR  +  FA 

or  specificity  =  1  -  a  (6) 

and  /?=  M  =\- sensitivity  =  power  of  the  test  (7) 
M  +  H 

with  sensitivity  =  \  -  f  (8) 

Note,  also  in  Figure  (3):  d'=— — —  (9) 

cr 

(3.  A)  Optimality  Measures  for  decision  making 

Three  types  of  optimality  in  decision  making  are  briefly 
discussed  herein.  The  first  two  types  of  decision  making  can 
be  developed  using  the  d’  and  c  values  in  Figure  (3). 

Optimality  Test  -1  -  Neyman-Pearson  (fix  a,  minimize  13) 

From  Figure  (3)  and  the  classical  maximum  likelihood 
ratio  test,  it  is  well  known  [2]  that  using  a  decision  rule  based 
on  the  likelihood  ratio  is  optimal  in  the  Neyman-Pearson 
sense.  This  means  for  a  fixed  a  (type  1  error),  the  minimum 
type  2  error  (/>  in  Figure  (3))  is  realized.  The  decision  rule 
employs  the  ratio  of  the  probability  density  functions  in  Figure 
(3).  The  decision  rule  is  to  select  choice  Hi  over  choice  H0  if 
(threshold  =  unity  in  Figure  (3)): 

X  =  (fj/fo)  >  Threshold  =  1  (10) 

The  second  and  well  accepted  test  of  optimality  involves 
the  area  under  the  ROC  (receiver  operator  characteristic) 
curve. 

Optimality  Test  2  -  Area  under  an  ROC  curve 

In  the  fields  of  medical  decision  making  and  diagnostisity 
in  general,  a  concept  called  “discriminability”  is  highly  touted 
in  these  fields.  Herein  a  brief  description  of  the  ROC  curve  is 
provided  to  be  instructive.  Figure  (4)  shows  the  ROC  curve  as 


(2.D)  -  Extant  Data  in  Human  Decision  Making. 

In  the  literature,  there  is  extensive  evidence  indicating  that 
human  decision  making  is  a  time-varying  process.  In  studies 
on  fatigue  or  vigilance,  this  effect  is  seen  quite  often.  For 
example  [3-5]  it  is  shown  that  d’  in  Figure  (3)  may  vary  with 
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Figure  (4)  -  The  ROC  (Receiver  Operator  Characteristic)  Curve 

a  plot  of  sensitivity  versus  1 -specificity.  The  area  under  the 
ROC  curve  (discriminability)  accounts  for  the  general  ability 
of  a  statistical  test  to  include  the  two  types  of  errors  (misses 
and  false  alarms).  It  represents  an  overall  measure  of  the 
efficacy  of  the  test.  This  is  precisely  what  PPP  and  NPP  can 
bring  to  the  decision  making  process. 

From  Figure  (4),  the  efficacy  of  the  test  (discriminability) 
is  the  total  area  under  the  ROC  curve  which  is  desired  to  be 
maximized.  It  should  be  pointed  out  that  ROC  curves  are  not 
arbitrary  functions  and  have  certain  restrictive  properties. 

The  next  optimality  test  will  generalize  signal  detection 
theory  to  dynamic  systems  which  is  appropriate  if  the 
moments  p(t)  and  a(t)  have  time  variations. 

Optimality  Test  3  -  The  Wald  Sequential  Test 

The  classical  Wald  sequential  algorithm  test  involves  two 
types  of  optimality  and  provides  an  entry  into  dynamic 
decision  making  [9].  Figure  (5)  portrays  a  dynamic  decision 
making  situation  in  which  the  ROC  curve  may  be  changing 
with  time. 


Hit  Rate  =  sensitivity  =  Sn  = 

in 


II +M 


(1.0,  1.0) 


False  Alarm  Rate  1.0 

FA 

-  1  -  specificity  - - 

FA+CR 


Figure  (5)  —  Adapting  ROC  curves  as  System  Parameters  Vary 

In  Figure  (5)  the  d’  and  c  parameters  may  vary  with  time 
and  the  goal  is  to  adjust  the  decision  rule,  accordingly.  To 
complete  this  background  information,  the  optimal  Wald 
sequential  decision  test  can  be  briefly  described  as  follows  [9]: 

A  statistical  measure  y(t)  (likelihood  ratio)  is  computed  in  real 
time  (Wald  test),  where 

A(t)  =  log^-  <  y(t)  <  log^=B(t)  (11) 
1  -a  a 


The  decision  rule  is:  If  the  determined  y(t)  >  B(t),  then  select 
hypothesis  H  in  Figure  (3).  However,  if  y(t)  <  A(t),  then 
choose  H0.  If  neither  these  conditions  are  true,  then  collect 
more  data. 

The  beauty  of  the  Wald  test  is  that  if  a  or  /?  vary  with  time, 
then  the  statistical  parameters  A(t)  and  B(t)  adjust  accordingly. 
The  decision  making  process  then  modifies.  It  is  shown  [9] 
that  the  Wald  test  still  maintains  two  types  of  optimality:  (1)  It 
is  still  Neyman-Pearson  optimal  in  the  sense  that  for  a  fixed  a 
(type  1  error)  then  P  (type  2  error)  is  minimized,  as  well  as  (2) 
the  algorithm  will  converge  in  minimum  time  to  a  decision. 
What  this  means  is  that  a  decision  is  reached  in  about  50%  of 
the  time  [9]  as  compared  to  a  static  maximum  likelihood 
decision  making  process  with  constant  values  of  a  and  />. 

Thus  dynamic  decision  making  with  time  varying  variables 
has  significant  advantages  over  static  methodologies. 

Finally,  the  concepts  of  PPP  and  NPP  can  now  be 
generalized  into  dynamic  decision  making  employing  the 
background  materials  developed  so  far. 


IV.  Adapting  Decision  Making  with  PPP  and  NPP 

Additional  facts  on  the  derivation  and  relationships  of  PPP 
and  NPP  to  dynamical  statistical  hypothesis  testing  are 
detailed  here.  A  brief  description  of  the  pertinent  steps  will  be 
presented  here.  Equations  (12-19)  now  show  key  relationships 
between  the  classical  quantities  CR,  M,  H,  and  FA  from  Figure 
(1)  to  the  PPP  and  NPP  parameters  as  well  as  sensitivity  (Sn) 
and  specificity  ( Sp ). 


CR=M\ 


NPP  ) 
1  -  NPP  J 


(12) 


( 1  -  NPP  \ 

M  =  CR 

1,  NPP  J 

(13) 

H  =  (FA)[  PPP  1 
\l-PPPJ 

(14) 

f\~ppp\ 

FA-{H\  PPP  ) 

(15) 

CR={FA)-Sp- 
1  -  Sp 

(16) 

FA-(CR)'-Sp 

Sp 

(17) 

H  =  (M )  Sn 

1  -Sn 

(18) 

M  =  (H) 


1  -Sn 
Sn 


(19) 
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The  following  key  result  has  its  details  contained  in  Appendix 
A  and  can  be  described  via  Theorem  1 : 

Theorem  1:  Sensitivity  ( Sn )  and  specificity  (Sp)  can  be 
related  to  PPP  and  NPP  through  the  following  formula: 


(  Sn  \ 

f  c  X 

Sp 

= 

f  ppp  ) 

f  NPP  } 

V  1  -S11  ; 

U-SpJ 

li  -PPP) 

[\-NPPJ 

Proof:  Appendix  A  derives  the  result  in  equation  (20). 

The  second  major  result  is  to  extend  the  Wald  sequential  test 
from  equation  (11)  into  a  form  for  the  use  of  PPP  and  NPP. 
With  this  dynamic  decision  making  criterion  the  Neyman- 
Pearson  property  of  optimality  is  still  preserved  and  the 
minimum  time  for  a  decision  to  be  made  is  still  left  intact. 
Theorem  2  presents  this  main  result: 

Theorem  2: 

The  equivalent  to  the  Wald  sequential  test  of  equation  (11)  can 
be  written  in  terms  of  PPP  and  NPP  via  the  following  method: 

The  optimal  Wald  sequential  algorithm  can  be  modified  to 
include  the  PPP  and  NPP  terms.  For  the  test  statistic  y(t) 
(likelihood  ratio)  in  equation  (11),  the  new  regions  for 
testing  now  include: 


log 


p\(\-PPPY\-NPP 


PPP 


log 


NPP 

PPP 


<_  y(t)  < 


NPP 


P 


\- all- PPP  \  \- NPP 


(21) 


Proof:  Appendix  B  outlines  the  derivation  of  equation  (21). 
Equation  (21)  is  equivalent  to  equation  (11)  and  still  preserves 
the  two  types  of  optimality  as  enjoyed  by  the  Wald  sequential 
test. 


Remarks:  The  formulation  (21)  has  new  advantages  over 
prior  methods.  It  is  known  that  the  parameters  d’  and  c  of 
Figure  (2)  may  vary  with  time.  From  equation  (21)  it  is  seen 
that  the  decision  parameters  a  and  p  are  known  apriori.  Also 
the  experimental  parameters  PPP  and  NPP  adapt  with  time 
since  the  d’  and  c  values  are  known  to  change  but  are  easily 
measured  in  real  time.  Thus  by  using  equation  (21)  rather  than 
(1 1)  provides  a  dynamic  decision  making  procedure  which  can 
be  modified  with  time  having  an  emphasis  on  the  time  varying 
PPP  and  NPP  variables.  The  result  in  equation  (21)  still 
enjoys  the  same  optimality  properties  as  in  the  Wald 
sequential  algorithm  in  equation  (11)  but  is  cast  within  the 
framework  of  the  PPP  and  NPP  variables. 


V.  SUMMARY  AND  CONCLUSIONS 


Sensitivity  and  specificity  using  a  PPP  and  NPP 
formulation  demonstrate  that  a  dynamic  Wald-type  sequential 
algorithm  can  be  synthesized  depending  only  on  the  PPP, 
NPP,  a  and  /?  values.  From  equation  (20),  several  new 
viewpoints  on  dynamic  decision  making  can  be  obtained.  For 
example,  if  PPP  may  be  constant  and  NPP  may  be  decreasing 
(cf.  [4]),  then  both  sensitivity  and  specificity  may  change. 


This  because  from  equations  (1,3)  PPP  and  Sn  are  different 
quantities.  The  new  viewpoint  allows  us  to  drill  down  on  the 
decision  making  when  overall  measures  like  (d’,  c,  Sn  and  Sp 
may  be  changing  with  time.  Also,  the  new  dynamic  decision 
making  algorithm  presented  herein  preserves  the  two  types  of 
statistical  optimality  provided  by  the  Wald  method  but  uses 
more  modern  procedures  such  as  PPP  and  NPP. 
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Appendix  A  -  Derivation  of  Equation  (20) 


Starting  with  the  equation  (4),  multiplying  through  by 
(CR+M)  yields: 


(  CR+M)  NPP  =  CR 

(A.  1) 

or 

CR  (l-NPP)  = 

-M  (NPP) 

(A-2) 

which  gives  rise  to  the  following  two  relationships: 

s; 

ii 

a 

(  NPP  ) 

(A.3) 

[l-NPP  j 

M  =  CR\ 

(l-NPP ) 

(A.4) 

(  NPP  ) 

Starting  with  equation  (3),  multiply  through  by  ( H+FA )  to 

yield: 

(H+FA)  PPP  =  H 

(A.5) 

or 

H (l -PPP)  =  FA  (PPP) 

(A.6) 
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which  arrives  at  the  following  two  relationships: 

PPP  \  (A. 7) 


H  =  {FA)  I 
FA  =  (H) 


1  -PPP 
l -PPP 
PPP 


CR  =  (FA) 
FA  =  (CR) 


1  -Sp 
l-Sp 
Sp 


(A.  11) 
(A.  12) 


Finally,  starting  with  equation  (1),  multiplying  through  by 
(H+M)  yields: 

Sn  (. H+M)  =  H 


or 


H(\-Sn)  =  Sn  (M) 


which  gives  rise  to  the  following  relationships: 
H=(M)~  Sn 


M  =  (Ft) 


1  —  Sn 
1  -Sn 
Sn 


To  now  develop  the  dependency  between  PPP,  NPP,  Sn  and 
Sp,  set  CR=CR  using  equations  (A. 3)  and  (A.l  1)  yielding: 
NPP 


M 


1  -NPP 


=  FA 


r  Sp  ^ 
1  -  Sp 


CR 


l -NPP 


=  H 


1  -Sn 
Sn 


V  NPP  J 

Following  the  same  procedure  by  setting  H=FI  in  equations 
(A. 7)  and  (A.  15)  which  yields: 


FA\ 


PPP 


=  M\  (A. 19) 

ki-ppp; 

Finally  setting  FA=FA  in  equations  (A. 8)  and  (A.12)  results 
in  the  following  relationship: 


1  -PPP  , 

H\  |  =  CR 

PPP 


1  -Sp 
V  SP 


The  relationships  (A.  17)-(A.20)  can  be  further  reduced. 
Computing  the  relationship  (CR)/H  from  both  equations 
(A.  18)  and  (A.20)  yields: 


1  -PPP 

) 

CR  _ 

l  Sn  J  l 

PPP 

J  (A.21) 

H  “ 

fl -Sp' 

[  NPP  J 

l  Sp  j 

fl-Sn') 

(\-PPP ^ 

fl-  NPP ") 

l  Sn  J 

l  Sp  J 

l  PPP  ) 

1  NPP  J 

To  check  why  this  is  true,  similar  calculations  can  be  made 
using  both  equations  (A.  17)  and  (A.  19)  for  (FA)/M: 


(A.23) 


f  S"  ) 

(  NPP  \ 

(A.8) 

FA 

yi-Sn) 

_U  -NPP) 

To  continue,  start  with  equation  (2)  and  multiply  through  by 
(< CR+FA )  to  yield: 

Sp  (CR  +  FA)  =  CR  (A.9) 

or:  CR  (  1  -Sp )  =  Sp  FA  (A.10) 

which  results  in  the  following  relationships: 

Sp 


M 


(  PPP  f 

(  c  \ 

Sp 

U  -PPP) 

Cross  multiplying  implies  the  reciprocal  of  equation  (A.22)  is 
true  which  further  validates  this  approach  because  if  things 
are  equal,  their  reciprocals  should  also  equate: 


(A. 24) 


(  Sn  ) 

f  C  ^ 

Sp 

f  PPP  ■) 

(  NPP  \ 

U  -  Sn  J 

{l-Spj 

[l  -PPP) 

u  -NPP) 

(A.13) 
(A.  14) 

(A. 15) 
(A.  16) 


Appendix  B  -  The  Wald  Algorithm  for  PPP  and  NPP 

Starting  with  the  Wald  test  of  equation  (11),  where  a  and 
/?  are  known  or  measured  (possibly  time  varying). 

J-  <  yO  <  iog^  (B.i) 


log 


1  -a 


a 


where,  from  equations  (6)  and  (8) 


and 


a  =  1  -  Sp 
P=\-Sn 


(B.2) 

(B.3) 


(A.  17) 


Thus  the  Wald  optimality  criterion  (B.I)  can  be  written: 

l-Sn  <  y(t)  <  j 
Sp  1  -Sp 


But  it  is  known  that 


Now  set  M=M  in  equations  (A. 4)  and  (A.  16)  which  results  in: 


(A. 18) 


f  Sn  f 

^  c  \ 

Sp 

(  PPP  \ 

(  NPP  \ 

U-ShJ 

ll -Sp) 

v  1  -  PPP  J 

U  -NPP) 

Thus 


f  1  -Sn) 

f  Sn  f 

fl-NPPf 

{  SP  J 

ll-^J 

{  PPP  j 

{  NPP  j 

f  1  -Sn') 

'i -/n 

(1  -PPP) 

fl-NPP) 

{  sP  J 

1  a  J 

V  PPP  J 

V  NPP  J 

(B.4) 

(B.5) 

(B.6) 

(B.7) 


(A.20) 


Also  starting  with  equation  (B.5)  yields: 


Sn  ^ 

f  PPP  f 

(  NPP  f 

l  Sp  J 

\\-ppp) 

[l-NPP  J 

(B.8) 


Which  can  be  written  as: 


(  Sn  'j 

r  fi  \ 

(  PPP  \ 

f  NPP  \ 

U  ~a) 

ll -PPP) 

ll -NPP) 

Cross  multiplying  implies  the  following  classical  dependence 
between  PPP,  NPP,  Sn  and  Sp: 


(A.22) 


(B.9) 


Flence  the  optimal  Wald  sequential  algorithm  can  be 
modified  to  include  the  PPP  and  NPP  terms  as  follows:  For 
the  test  statistic  y(t)  assuming  constant  a  and  ft  values  with 
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NPP  and  PPP  possibly  changing  with  time,  the  new  regions 
for  testing  now  become: 


log 


i  -ppp'] 

fl-NPP\ 

[l « 1 

ppp  j 

t  NPP  J 

<_log 


p  Y  ppp 


1  -  a  A  1  -  PPP  A  1  -  NPP 


<  y(t) 

NPP 


(BIO) 
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Abstract 

New  applications  of  approximate  entropy  are  currently  ongoing  and  being  analyzed  with  innovative 
and  novel  data  being  generated  within  AFRL.  This  paper  reviews  some  extant  studies  in  the  area  of 
Human  Machine  studies  and  describes  the  basics  of  this  powerful  use  of  an  infonnation  theoretic 
measure.  The  goal  is  to  detect  if  the  state  of  a  human  may  have  been  compromised  by  fatigue,  or  other 

stressors,  utilizing  noninvasive  sensory  measurements. 


I.  Introduction 

New  investigations  at  the  711  HPW  WPAFB  have  been  focused  on  the  use  of  a  measure  of 
uncertainty  in  data  termed  “approximate  entropy.”  This  new  interest  in  research  in  the  USAF  involves 
studies  of  humans  under  fatigue  induced  by  loss  of  sleep  [1]  or  other  stressors  including  emotional  stress 
[2],  Herein  will  be  presented  some  of  the  basics  of  the  concepts  underlying  the  calculation  of  approximate 
entropy.  Prior  work  in  this  area  [3]  has  considered  detennining  when  loss  of  consciousness  would  occur 
during  high  acceleration  stress.  The  study  in  [3]  was  based  on  an  exogenous  physical  stress  (G  or 
acceleration  forces)  on  a  human.  In  [4],  however,  a  mental  trauma  involving  cognitive  stress  (a  high 
workload  task)  was  also  investigated  using  the  approximate  entropy  measure.  Periods  of  high  and  low 
task  difficulty  could  be  gleaned  from  the  performance  of  the  subjects  and  the  approximate  entropy 
measure. 

II.  Principles  of  Approximate  Entropy 

The  concepts  involving  approximate  entropy  are  simple.  In  Figure  (1)  a  signal  s(t)  =  Si(t)  is 
compared  to  itself  displaced  one  sample  s(t+At)  =  S2(t).  The  entropy  is  the  normalized  amount  of  disorder 
between  the  two  adjacent  signals  Si(t)  versus  S2(t)  or  s(t)  versus  s(t+At).  If  the  approximate  entropy  is  low, 
this  disorder  is  minimal.  If  this  disorder  or  variation  is  high,  then  the  normalized  entropy  may  reach 
values  of  1.0  or  larger.  In  references  [3]  and  [4]  more  technical  details  are  included.  Appendix  A  includes 
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reference  [3]  involving  physical  stress  on  humans  for  informative  purposes.  Appendix  B  includes 
reference  [4]  with  a  cognitive  stress  investigation. 


How  is  Approximate  Entropy  Calculated? 


(1)  Correlate  a  signal  Si(t)  with  itself  S2(t) 

(2)  Slide  window  along  signal’s  data  points 

(3)  Analyze  signal  entropy  (disorder) 

=-2>log(P) 

Figure  (1)  -  Determination  of  the  Approximate  Entropy  Measure 
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Appendix  A  -  Reference  [3]  included  for  information  purposes 


APPROXIMATE  ENTROPY-AN  ASSESSMENT  TOOL  FOR  SYSTEM  COMPLEXITY  AND 

UNCERTAINTY 

D.  W.  Repperger1,  W.  B.  Albery1,  L.  D.  Tripp2 


1  Air  Force  Research  Laboratory,  AFRL,  WPAFB,  Ohio  45433,  USA 
2  General  Dynamics,  Dayton,  Ohio,  45433,  USA 


Abstract:  For  complex  systems  and  interactions,  a  means  of  determining  the  quality  of 
information  in  a  time  series  involving  data  generated  from  a  human-machine  system  is 
investigated.  This  recent  numerical  analysis  method,  widely  accepted  in  the  medical 
field,  termed  “Approximate  Entropy”  can  provide  a  compelling  means  of  uncovering 
irregularity  in  data  of  all  types.  Changes  in  the  state  of  a  human-machine  system  can  be 
quickly  gleaned  in  real  time,  on  line.  Application  of  this  measure  is  investigated  on 
tracking  performance  data  when  the  human-machine  system  is  known  to  be 
compromised  in  a  performance  sense. 

Keywords:  Computational  methods,  classifiers,  entropy,  discriminations,  dynamic  behavior. 

1.  INTRODUCTION 

Many  complex  Human-Machine  systems  may  have  a  sudden  change  in  state  and  it  is  desired  to 
detect,  in  real  time,  whether  a  measured  or  derived  signal  may  now  differ  from  its  prior  situation. 
Traditional  statistical  measures  (mean,  standard  deviation  or  other  moments)  cannot  always  capture 
the  change  in  some  signals  because  of  the  complexity  of  the  underlying  system  dynamics.  Pioneered 
by  Pincus  (1991-2),  discussed  in  Pincus  and  Kalman  (1997)  and  in  numerous  publications,  an 
interesting  new  measure  has  been  developed  to  help  detect  changes  in  the  “regularity”  of  a  time  series 
termed  “Approximate  Entropy”  (ApEn).  By  comparing  a  time  series  with  itself,  over  time,  this 
measure  of  system  irregularity  has  found  applicability  in  the  medical  field  for  discerning  differences 
in  heart  beat  either  due  to  sleep  cycle  or  disease,  Yeragani,  et  al.  (1998)  and  Pincus  (1992).  This 
procedure  has  applicability  in  distinguishing  honnonal  changes,  Gevers,  et  al.  (1998),  mood  swings, 
Pincus  (2003),  and  in  other  complex  applications  involving  very  intricate  signals.  Other  applications 
include  EMG  and  tremor  distinction  (Morrison  and  Newell,  2000),  EEGs  (Bruhn  et  al.,  2000), 
recognizing  epileptic  activity  (Diambra,  et  al.  1999)  and  for  discerning  cocaine  addiction  (Newlin  et 
al.,  2000).  Generalizing  this  concept  even  further,  Sugihar  (1990),  it  is  now  possible  to  determine  if 
the  underlying  complexity  of  a  system  state  may  have  changed  when  viewed  within  the  context  of  a 
control  system  by  observing  some  key  output  data.  This  methodology  has  particular  interest  in 
identifying  if  the  underlying  dynamics  of  a  process  may  or  may  not  exhibit  chaotic  behaviour  and  to 
predict  this  possibility. 

Approximate  entropy  is  based  on  the  simple  principle  that  if  a  time  series  signal  can  be  compared  to 
itself  (heart  beat  data,  for  example)  and  the  amount  of  the  disorder  in  a  comparison  or  change 
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(entropy)  between  relative  time  shifts  of  these  data  is  increasing,  this  is  probably  an  indication  of 
some  change  in  state.  This  differs  from  traditional  correlation  measures,  since  they  are  not  based  on 
information  theory  concepts  and  may  require  fixed  implicit  models.  The  use  of  ApEn  is  model 
independent  and  only  depends  on  the  real  time  data  series.  There  are  at  least  four  reasons  why 
approximate  entropy  provides  new  information  about  system  complexity  not  normally  derived  from 
typical  statistical  measures  (first  and  second  order  moments): 

(1)  If  data  are  noisy,  the  approximate  entropy  measure  can  be  compared  to  the  noise  level  in  the  data 
to  determine  what  quality  of  true  information  may  be  present  in  the  data. 

(2)  If  the  data  have  an  artifact,  this  does  not  impact  the  approximate  entropy  measure  as  much  as  it 
would  affect  typical  first  and  second  order  statistical  moments  from  the  data. 

(3)  Approximate  entropy  can  be  designed  to  work  for  small  data  samples  (n  <  50  points)  and  can  be 
applied  in  real  time,  on  line.  Thus  changes  in  the  state  of  a  physical  process  may  be  quickly 
detennined. 

(4)  For  pure  stochastic  processes,  approximate  entropy  will  become  practically  infinite. 

Thus,  the  quality  of  the  information  in  a  signal  can  then  be  quantitatively  evaluated  by  comparing  the 
entropy  level  of  the  measured  signal  with  its  underlying  (non  random)  signal  component. 


To  test  the  concept  of  approximate  entropy  on  perfonnance  data,  an  extensive  database  exists 
involving  pilots  performing  cognitive  tasks  while  simultaneously  being  subjected  to  high  acceleration 
stress,  Tripp  (2001).  The  volunteer  subjects  were  required  to  perfonn  both  a  math  computation  task 
as  well  as  manual  tracking  in  a  position-control  mode.  Data  were  recorded  prior  to  the  event  of  loss  of 
consciousness,  during  the  period  the  subjects  went  unconscious,  and  afterwards  in  the  post  recovery 
period.  The  approximate  entropy  measure  was  calculated  on  the  tracking  data  to  quantify,  objectively, 
the  amount  of  disorder  in  the  subject’s  response  during  different  time  epochs  of  the  experiment. 

2.  EXPERIMENTAL  SCENARIO 

At  the  Air  Force  Research  Laboratory,  WPAFB,  Ohio,  the  data  presented  here  partially  involved  a 
three  axis  motion  simulator  that  is  used  to  detennine  a  pilot’s  response  to  acceleration  stress.  The 
system  shown  in  Figure  1  has  a  5.8  meter  radius  with  a  large  spherical  cab  and  can  create  a  force  of 
20  G  at  a  rotational  velocity  of  56  RPM.  Such  a  system  weights  163,000  kilograms. 
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Figure  1  -  The  Dynamic  Environment  Simulator  to  Produce  Acceleration  Stress 

3.  OBJECTIVE 

The  objective  is  to  test  if  the  measure  of  uncertainty  or  irregularity  known  as  approximate  entropy  can 
be  a  valuable  tool  in  assessing  whether  a  pilot  is  compromised  in  a  cognitive  sense.  The  underlying 
assumption  is  that  under  high  acceleration  stress  conditions,  human  data  are  extremely  more  variable 
and  thus  should  show  higher  entropy.  In  addition,  if  the  time  rate  of  increase  of  entropy  is  positive 
going  into  the  loss  of  consciousness  event  and  negative  coming  out  of  the  loss  of  consciousness  event, 
then  the  rate  of  change  of  approximate  entropy  may  be  a  valuable  prediction  tool  to  ascertain  shifting 
cognitive  state.  Also,  if  the  rate  of  change  of  approximate  entropy  is  decreasing,  this  may  be  a 
valuable  indicator  that  the  mission  capability  of  the  pilot  has  returned  to  a  normal  level  and  the  effect 
of  the  untoward  event  (unconsciousness)  may  no  longer  impact  the  mission  effectiveness  of  the  pilot. 

4.  HYPOTHESIS 

The  null  hypothesis  that  we  wish  to  reject  is  that  the  approximate  entropy  metrics  (magnitude  and 
time  rate  of  change  of  ApEn)  will  not  significantly  vary  during  known  changes  of  the  cognitive  state 
of  the  pilot  being  stressed.  Only  the  performance  tracking  data  will  be  used  to  ascertain  the  cognitive 
state  of  the  pilot. 


5.  METHODS 

Data  from  16  USAF  pilots/subjects  were  collected  as  they  performed  a  compensatory  tracking  task 
and  carried  out  a  math  computation  exercise  in  a  large  centrifuge  motion  simulator.  These  subjects 
perfonned  these  tasks  under  high  acceleration  stress.  The  acceleration  stress  increased  until  the 
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subjects  became  unconscious  (GLOC  or  G  loss  of  consciousness).  There  were  3  or  4  data  days  for 
each  subject. 


6.  APPARATUS 

Two  centrifuges  were  used  in  this  study.  The  DES  centrifuge  at  WPAFB  was  described  in  section  2 
of  this  paper.  A  second  centrifuge  at  Brooks  Air  Force  Base,  Texas  was  employed  to  generate  similar 
acceleration  profiles.  The  same  performance  tasks  (math  computation  and  manual  tracking)  were  used 
with  both  motion  simulators  on  qualified  human  subjects. 

7.  EXPERIMENTAL  DESIGN 

Both  genders  of  healthy  USAF  pilots/subjects  participated  in  this  experiment.  The  end  point  condition 
was  the  loss  of  consciousness.  A  recovery  period  followed  the  GLOC  event. 

8.  RESULTS 

Figure  2  shows  the  stressor  (lower  plot  -  G  acceleration  level  versus  time)  and  the  root  mean  square 
tracking  error  (upper  plot  in  Figure  2)  for  one  subject  on  the  first  of  his  four  days  when  data  were 
collected.  The  events  of  pre-GLOC  (prior  to  G  loss  of  consciousness),  the  period  of  complete 
incapacitation,  and  the  post  recovery  period  are  indicated.  Figures  3a-b  is  a  similar  plot  of  the 
approximate  entropy  function  versus  time  with  some  of  these  same  event  markers  noted  based  on  10 
second  intervals  of  the  data.  One  observes  in  Figure  2,  from  the  performance  data  (top  plot  of  Figure 
2),  that  prior  to  the  GLOC  event,  a  pilot’s  behaviour  was  manifested  by  a  high  variation  in  the 
tracking  error  signal.  Coming  out  of  the  GLOC  event,  in  the  same  diagram,  high  variability  in  the 
tracking  error  also  exists  for  40  or  more  additional  seconds.  Eventually  the  data  show  a  reduction  in 
the  variability  (t  >  100  seconds)  as  the  pilot  returns  to  a  cognitive  state,  typical  of  normal  tracking. 
Hence,  disorder  (high  levels  of  approximate  entropy)  in  the  tracking  error  signal  seems  related  to  a 
compromised  cognitive  state  of  the  pilot. 


The  GLOC  Event 


Figure  2  -  ERMS  and  G-Stress  Versus  Time 


Figure  3a  -  Subject  1  -  Day  1  -  3D  Surface  Plot  of  Approximate  Entropy 
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The  appendix  describes,  in  a  succinct  manner,  important  technical  aspects  of  this  measure,  approximate 
entropy,  as  compared  to  alternative  metrics  involving  time  series  data.  With  reference  to  Figures  3a-b,  the 
approximate  entropy  ApEn  is  plotted  versus  time  on  the  x  axis.  One  thing  is  clear,  i.e.  prior  to  the  GLOC 
event  and  during  the  post  GLOC  event,  that  the  absolute  value  of  ApEn  and  the  magnitude  of  the  rate  of 
change  of  ApEn  with  respect  to  time  are  extremely  high.  In  Figures  3a-b,  the  almost  zero  value  of  ApEn 
near  the  fifth  time  epoch  (t  =  50  seconds)  needs  to  be  explained.  Referring  to  Figure  2,  during  the  time 
period  (t  s  [25,35])  seconds,  the  RMS  tracking  error  is  constant  (worst  case)  because  the  subject  is 
unconscious.  There  is  little  variation  in  his  relative  behavior  (as  measured  by  tracking  performance, 
which  is  not  changing),  hence  ApEn  — >  0  in  this  interval.  Therefore  to  use  a  measure  such  as  ApEn  to 
predict  or  indicate  the  likelihood  of  the  pilot  to  become  unconscious,  one  must  take  into  account  the 
magnitude  of  ApEn  as  well  as  it’s  time  rate  of  change.  Figure  4  depicts  a  majority  voting  classifier 
system  that  could  use  this  information  to  predict  the  cognitive  state  of  the  pilot  using  as  inputs  ApEn 
measures  and  tracking  performance.  Since  brevity  must  be  the  style  in  this  paper,  we  only  report  on  the 
data  analysis  measure  of  ApEn  with  respect  to  the  known  cognitive  state  of  the  pilot. 
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Figure  4  -  A  Majority  Voting  Scheme  to  Predict  the  Pilot’s  Cognitive  State 


Table  1  illustrates  the  mean  (|u)  and  standard  deviation  (SD)  values  of  ApEn  and  its  time  derivative 
fifteen  seconds  prior  to  the  GLOC  event.  Also  in  Table  1  are  these  respective  values  100  seconds  post 
GLOC  for  comparison  purposes.  The  data  are  averaged  across  16  subjects  with  3-4  data  days  for  each 
subject.  Figure  5  illustrates  these  mean  and  SD  values  averaged  over  the  16  subjects  and  for  those  data 
days  with  pre  GLOC  data  at  15  or  more  seconds  prior  to  the  incapacitation  event. 

Table  1  -ApEn  and  its  time  derivative  15  seconds  Pre  GLOC  and  100  seconds  Post  GLOC 

Pre  GLOC  Post  GLOC 

ApEn  (d/dt)  ApEn  ApEn  (d/dt)  ApEn 
p  SD  p  SD  p  SD  p  SD 

.43  .07  .013  .006  .32  .13  .01  .006 

Table  2,  displays  the  results  of  a  one-way  ANOVA,  one  factor  for  Magnitude  of  ApEn  and  Table  3 
displays  these  values  for  the  magnitude  of  the  time  rate  of  change  of  ApEn.  Please  note  the  time  samples 
are  one  second  apart  for  the  data  displayed  in  Figure  2.  Figure  6  shows  the  results  of  the  JMP  (developed 
by  the  SAS  Institute)  analysis  for  comparison  of  means  (Tukey-Kramer  test). 
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Magnitude  of  ApEn  prior  and  postGLOC 


ApEn  Magnitude 
Values  Normalized 


1 


2 


Series  1  -  Mean  and  SD  Values  Pre  GLOC 
Series  2  -  Mean  and  SD  Values  Post  GLOC 


Time  rate  of  change  of  ApEn  for  Pre  versus  post  GLOC  -  Means  and  SD 


Figure  5a  -  Magnitude  of  ApEn  for  Pre  versus  Post  GLOC 


Series  2  are  the  Mean  and  SD  values  Post  GLOC 

Figure  5b  -  Magnitude  of  d/dt(ApEn)  for  Pre  versus  Post  GLOC 


Table  2  -  One  Way  ANOVA  for  Mag,  of  ApEm 
Source  DoF  £  Squares  Mean  Squares  F  Ratio  Prob>F 
GLOC  1  0.2887  0.2887  20.7408  <.0001 

Error  70  0.9742  0.0139 

C.  Total  71  1.2629 
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Table  3  -  One  Way  ANOVA  for  Mag  of  d/dt  (ApEn) 

Source  DoF  X  Squares  Mean  Squares  F  Ratio  Prob>F 


GLOC 


1  0.00024 

Error  70 
C.  Total 


0.00024  6.3978 

0.0026  0.000037 

71  0.00284 


.0137 


Figure  6  -  JMP  Analysis  for  d/dt  (ApEn)  Data  -  Pre  and  Post  GLOC. 

9.  DISCUSSION 

Preliminary  results  indicate  that  prior  to  GLOC  the  approximate  entropy  function  is  both  high  and 
shows  increasing  positive  rates  with  respect  to  time.  In  the  recovery  period,  the  approximate  entropy 
function  seems  to  show  decreasing  rates  as  the  pilot’s  variability  in  his  performance  seems  to  be  more 
reduced.  The  overall  results  of  this  study  across  all  16  subjects  with  three  to  four  days  of  testing  for 
each  subject  are  reported. 


10.  CONCLUSION 

A  study  of  the  efficacy  of  using  approximate  entropy  to  predict  the  cognitive  state  of  the  pilot  is 
conducted  on  human  subjects  known  to  be  compromised  in  a  performance  sense  due  to  loss  of 
consciousness.  Studies  are  ongoing  to  predict  the  loss  of  consciousness  prior  to  the  event  by 
examining  the  rate  of  change  of  approximate  entropy  in  real  time,  on  line,  using  other  signals  such  as 
physiologically  based  data  including  oxygen  perfusion  to  the  subject’s  brain,  etc. 
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APPENDIX  from  [3]  -  TECHICAL  DETAILS  REGARDING  APPROXIMATE  ENTROPY 

Approximate  entropy  is  sometimes  termed  a  “regularity  measure”  which  quantifies  the 
unpredictability  of  fluctuations  in  a  time  series,  e.g.  in  instantaneous  heart  rate  signal.  High  levels  of 
ApEn  (highly  irregular)  reflect  the  likelihood  that  “similar”  patterns  of  observations  will  NOT  be 
followed  by  additional  “similar”  observations.  Thus  a  more  complex  process  has  high  levels  of  ApEn 
and  small  values  of  ApEn  imply  predictable  (repetitive)  patterns  are  inherent  in  the  data.  Uncertainty 
(or  high  system  complexity)  is  related  to  high  levels  of  ApEn. 

Conversely,  low  values  of  ApEn  indicate  predictability  of  a  time  series.  Given  a  series  of  s(t) 
measurements  (for  this  case  n  <  50)  s(  1),  s(2),  . . .  s(n),  equally  spaced  in  time,  the  ApEn  of  this  data 
series  depends  on  two  key  parameters:  m  and  r.  m  is  an  integer  that  represents  the  length  of 
compared  runs  (a  window  or  how  many  data  samples  the  two  series  differ)  and  r,  effectively, 
represents  a  filter.  Typically  m  =  1  or  2  which  distinguishes  the  two  series  and  r  is  a  tolerance 
measure  (criterion  of  similarity).  The  next  steps  provide  the  procedure  to  compute  the  ApEn  versus 
time: 
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Step  1:  Form  a  sequence  of  vectors  x(l),  x(2), x(n),  where  each  x(i)=[s(i),  s(i+l),  ...,s(i+m-l)]. 


Step  2:  Use  the  sequence  x(l),  x(2),  . . x(n)  to  construct,  for  each  i,  l<i<n-m+l,  Qm(  r  )  =  (the 
number  of  x(j)  such  that  d[x(i),  x(j)]  <  r)  /  (n-m-1).  The  distance  metric  d  satisfies: 

d[x(i),x(j)]  =  max  |s(i+k-l)  -  s(j+k-l)|  (A.l) 

for  k=l,2,...,m 


Hence  d  represents  the  distance  between  the  vectors  x(i)  and  x(j),  given  by  the  maximum  of  their 
respective  scalar  components.  The  next  step  defines  the  logarimithic  entropy  measure: 


Step  3: 

4>m(  r  )  =  (n-m+1)'1  I  In  C,m(r) 

resulting  in 

ApEn  =  c|)m(  r  )  -  c^m+1(  r  ) 


n-m+1 

i=l 


(A. 2) 

(A. 3) 


We  now  show  a  relationship  of  ApEn  to  some  of  the  objective  means  of  quantifying  chaos  in 
nonlinear  dynamical  systems.  Three  popular  means  of  quantifying  chaos  are  briefly  discussed  here 
which  include  (1)  Lyapunov  Exponents,  (2)  Embedding  Dimensions,  and  (3)  K-S  Entropy.  Large 
values  of  ApEn  can  be  related,  e.g.,  to  a  positive  Lyapunov  exponent. 

(1)  Lyapunov  Exponents: 

In  a  chaotic  system,  the  Lyapunov  Exponent  represents  the  exponential  rate  of  separation  of  adjacent 
trajectories.  Consider  a  nonlinear  system  near  a  fixed  point  (xo)  where: 

d/dt  x(t)  =  f(x(t))  (A.4) 

Taking  the  Taylor  series  expansion  near  xo  results  in: 

d/dt  x(t)=  0  =  f(x)  =  f(xo)  +  (x-xo)  df/dx  +...  (A. 5) 

Define  a  new  variable  z=x-xo.  Then  proximal  to  the  fixed  point:  d/dt  z(t)  =  z  [df/dx]  |xo  which  has 
solution:  z(t)  =  z(0)  e  ,  where  the  Lyapunov  exponent  X  =  [df/dx]  |x0  is  a  characteristic  value  of  the 
fixed  point.  Thus  two  adjacent  trajectories  are  attracted  to  each  other  if  they  approach  the  fixed  point 
and  X  <0.  However,  if  two  adjacent  trajectories  repel  each  other,  then  X  >  0  and  this  is  typical  of 
chaotic  behaviour.  Thus  if  the  average  Lyapunov  exponent  is  positive  near  a  fixed  point,  this  is  a 
viable  definition  of  chaotic  behavior  since  nearby  trajectories  have  exponential  divergence  in  phase 
space  (Hilborn,  1994). 

(2) Embedding  Dimensions:  The  embedding  dimension  is  the  number  of  points  necessary  to  predict 
the  next  point  in  the  time  series  and  is  an  objective  measure  of  system  complexity.  Thus  highly 
complex  systems  would  have  a  high  measure  of  embedding  dimension  and  this  may  indicate  the 
potential  for  chaos  (Williams,  1997). 

For  K-S  Entropy:  The  Kolomogorov-Sinai  entropy  (Kolmogorov,  1958)  deals  with  the  probability 
that  a  given  trajectory  point  falls  within  some  particular  region  of  state  space.  Thus  a  chaotic  system 
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would  have  a  high  entropy  measure  since  it  may  have  a  nonzero  high  level  of  probability  of  accessing 
different  areas  of  the  state  space.  Pincus  (1991)  has  shown  the  distinction  between  ApEn  and  K-S 
entropy  by  noting  that  the  K-S  entropy  measure  (a  theoretical  construct)  can  be  defined  via: 

For  K-S  Entropy: 

lim  lim  lim  [(j)m(  r  )  -  cj)m+1(  r )]  =  At  h(p)  (A. 6) 

r— >0  m— >co  n— >oo 

where  h(p)  is  the  K-S  entropy.  Choosing  At=l  shows  the  relationship  and  distinction  between  ApEn 
and  the  K-S  entropy  measure.  In  other  words,  the  K-S  entropy  is  a  theoretical  construct  which 
requires  r  to  be  near  zero,  and  m  and  n  to  be  large  in  value.  ApEn,  on  the  other  hand,  may  have  r 
nonzero  (but  small)  and  finite  values  of  m  and  n.  Thus  ApEn  is  more  applicable  to  real  world 
situations  and  data  that  are  collected  in  an  experimental  scenario.  Other  methods  also  exist  for 
quantifying  chaos  include  fractal  dimension  and  correlation  dimension  (Hilborn,  1994). 

Finally,  we  describe,  and  show  by  example,  how  the  calculation  of  ApEn  can  be  conducted  online  in 
real  time  for  data  consisting  of  finite  samples.  To  physically  understand  the  key  parameters  r  and  m, 
Figures  7a-b  displays  a  time  series,  where  s(t)  is  the  signal  of  interest.  The  signal  s(t)  has  two 
constituent  signals  Si(t)  and  S2(t)  derived  by  shifting  the  data  one  data  sample  (m=l)  over  a  12  sample 
time  epoch  (n=12).  S2(t)  is  just  one  data  sample  (m=l)  shifted  to  the  right  from  Si(t).  Hence,  the 
tolerance  variable  r  would  be  gleaned  from  the  difference  between  Si(t)  and  S2(t)  as  indicated  in  Figure 
7b  (plotted  on  different  and  exaggerated  vertical  scales  for  illustrative  purposes).  Thus  the 
approximate  entropy  value  depends  heavily  on  the  n,  m,  and  r  values  as  well  as  the  regularity  of  the 
signal  s(t)  over  the  interval  in  which  it  is  being  examined. 

In  conclusion,  a  relationship  between  the  Lyapunov  exponent  and  the  prior  procedure  using  ApEn 
shows  a  further  hidden  connection  of  the  proposed  method  to  nonlinear  dynamics.  One  can  derive  the 
Lyapunov  exponent  from  real  data  (Wolf,  et  ah,  1985)  and  it  is  desired  to  distinguish  the  situation  of 
chaos  from  possible  measurement  error  in  a  time  series  (Sugihar  and  May,  1990)  using  real  time 
measurements.  To  show  how  to  employ  the  Lyapunov  exponent  method  with  real  data,  the 
assumption  is  first  made  that  the  data  are  sampled  equally  in  time.  The  data  samples  s(to), 
s(ti),s(t2),...,s(tn)  are  labelled  So,  si,  S2,  ...,  sn.  If  x  is  the  constant  time 

interval  between  samples,  then  for  some  integer  n,  the  following  relationships  exists: 

tn— 10  =  n  x  (A. 7) 

Now  a  system  will  behave  chaotically  if  the  divergence  of  nearby  trajectories  are  showing  exponential 
changes  in  their  differences,  i.e. 


* 

1 

II 

o 

T3 

(A.  8) 

di=  |  Xj+i-Xi+i| 

(A. 9) 

d2  =  |  Xj+2  -  Xi+2| 

(A.  10) 
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*• 

t  =  Time  in  seconds 
Figure  7a  -  Original  Signal  s(t)  to  be  analyzed  versus  time. 


Figure  7b  -  Constituent  signals  s,(t)  and  s2(t)  derived  from  s(t). 


is  assumed  to  exponentially  increase  (on  the  average)  as  n  gets  larger.  The  assumption  is  really  that: 

d„  =  d0e^n  (A.  12) 

and  the  quantity  X  can  be  determined  via 

X  =  (1/n)  In  (dn/d0)  (A.  13) 

and  if  X  is  positive,  then  the  behavior  would  be  chaotic.  A  graphical  means  to  determine  this  effect  is 
to  make  a  log  plot  of  the  respective  difference  ratios  (dn/do)  and  detennine  if  the  line  that  best  fits  the 
data  has  a  nonzero  slope. 

If  the  slope  of  the  line  is  nonzero  (and  statistically  different  from  zero)  then  X  ^  0  and  the  data  either 
represent  a  chaotic  system  (/>())  or  an  attractor  (a<0),  the  latter  of  which  would  indicate  convergence 
to  a  fixed  point.  Hence  X  is  similar  to  the  Lyapunov  exponent. 
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Appendix  B  -  Reference  [4]  included  for  information  purposes 


A  Real-Time  Measure  to  Study  Dynamic  Interactions  with  a  Visual  Display 
D.  W.  Repperger  and  J.  J.  Skelly 

Abstract:  As  operators  interrelate  with  complex  systems,  such  as  those  that  occur 
in  flying  or  operating  aircraft,  the  true  perfonnance  and  interaction  is  of  a  dynamic 
nature.  However,  most  analyses  of  human  interface  systems  employ  the 
assumption  that  the  display  is  static  and  does  not  change  with  time.  This  paper 
will  consider  a  real  time  measure  (approximate  entropy)  to  study  and  evaluate 
human-machine  interaction  as  events  in  the  environment  change  in  a  dynamic 
sense.  The  experimental  paradigm  to  evaluate  the  efficacy  of  the  proposed  real 
time  measure  will  be  borrowed  from  studies  investigating  the  effects  of  spatio- 
temporal  structure  on  cognition. 

Keywords:  Dynamic  behavior,  entropy,  man-machine  interfaces,  nonstationary. 

1.  INTRODUCTION 

In  many  physical  systems,  the  presumption  that  a  static  environment  captures  the  essence  of  the 
performance  of  the  human-machine  interaction  is  extremely  limiting.  This  is  equivalent  to  the 
assumption  of  linearity.  In  linear  systems,  the  transient  response  has  commonalty  with  the  steady 
state  response  through  key  parameters,  such  as  eigenvalues.  However,  for  nonlinear  systems,  this  is 
not  the  case  and  most  human  interactions  with  complex  systems  are  more  typically  nonlinear.  In  a 
companion  paper,  Repperger,  et  al.  (2004),  a  real-time  measure  of  entropy  or  disorder  has  been 
employed  to  analyze  data  from  an  experiment  when  it  was  known  that  pilots  became  unconscious  and 
their  behaviour  changed  dramatically.  The  real  time  measure  of  approximate  entropy  (ApEn  -  Pincus 
and  Viscarello  (1992),  Gevers  et  al.  (1998),  and  Pincus  (2003))  has  been  embraced  by  the  medical  and 
psychological  communities  to  show  that  certain  changes  in  real  time  data  can  indicate  an  alteration  in 
the  medical  or  psychological  state  of  a  patient.  This  is  not  unlike  humans  interacting  with  dynamical 
systems,  when  they  become  overloaded  and  task  perfonnance  changes,  accordingly.  This  study  will 
consider  the  real  time  measurement  of  irregularity  in  data  (approximate  entropy)  as  applied  to  key 
performance  parameters  from  a  spatial-temporal  investigation  on  human  perception.  It  is  known  that, 
in  scenarios  of  the  data  to  be  presented,  human  performance  can  be  compromised  by  the  procedure  at 
which  information  is  displayed  to  the  operator.  This  is  true  for  two,  almost  identical  performance 
tasks  to  be  offered,  where  the  complexity  of  the  task  (bits)  and  rate  of  presentation  (bits/second)  are 
identical  for  both  tasks,  but  their  dynamic  attending  characteristics  differ.  From  a  Fitts’  Law 
perspective  (Fitts,  1954)  the  two  tasks  provide  identical  difficulty  but  it  will  be  shown  in  the  sequel 
that  the  performance  of  the  subjects  is  markedly  dissimilar.  The  only  distinction  between  the  tasks  is 
that  they  require  different  dynamic  attention  resources  from  the  human  subjects.  This  experimental 
platform  is  ideal  for  evaluation  of  the  ApEn  metric  which  may  be  able  to  discern  differences  in 
operator  responses,  when  Fitts’  law  would  not  be  able  to  make  such  a  fine  distinction. 
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A  powerful  paradigm  to  investigate  the  effects  of  dynamic  information  structure  on  human 
perfonnance  can  be  seen  in  the  works  of  Skelly  and  colleagues  (Skelly  (2003),  Skelly  et  al.  (2000), 
and  Jones  and  Skelly  (1993)).  To  briefly  summarize  the  results  of  some  of  their  studies  relevant  to  this 
investigation,  the  environment  of  an  operator  is  made  to  change  both  in  a  temporal  and/or  spatial 
sense  with  respect  to  a  visual  display.  The  timing  and  spatial  display  of  this  information  to  the 
operator  may  affect  a  viewer’s  ability  to  allocate  attention  among  items  within  a  single  stream  of 
information  events  or  among  multiple  infonnation  streams.  Hence  these  spatio-temporal  structural 
properties  may  serve  to  facilitate  or  interfere  with  attending  to  the  specific  arrival  time  or  location 
associated  with  dynamic  visual  stimuli.  The  underlying  thesis  is  that  goal  directed  targeting  of 
attention  may  not  be  entirely  voluntary,  i.e.  attentional  targeting  may  be  influenced  by  joint  structural 
properties  associated  with  dynamic  stimuli,  and  thus,  attention  may  be  involuntarily  controlled,  at 
least  in  part  by  the  design  of  the  display  dynamics.  This  dynamic  attending  approach  is  a  biological 
view  on  how  we  allocate  our  attentional  resources  to  dynamic  information  in  the  environment.  There 
are  three  basic  assumptions  in  this  dynamic  approach  to  understanding  how  humans  interact  with 
complex  systems:  (1)  Attention  is  controlled,  in  part,  by  the  combined  spatial-timing  structure  of 
dynamic  visual  infonnation,  (2)  Attentional  energies  are  stimulated  and  “synchronized”  with  certain 
invariantly  occurring  dynamic  space-time  structures,  and  (3)  The  focus  of  the  attention  may  speed  up 
or  slow  down  when  there  are  abrupt  accelerations  or  decelerations.  It  has  been  clearly  demonstrated 
that  certain  spatial-time  patterns  are  productive  for  the  transfer  of  infonnation  from  the  display  to  the 
operator  and  other  patterns  can  be  very  counterproductive  for  this  task. 

2.  OBJECTIVE 


The  objective  in  this  investigation  is  use  the  experimental  platform  involving  dynamic  attending  with 
the  approximate  entropy  measure  on  performance  data  involving  a  discrete  sequence  of  tasks.  A 
causality  will  be  shown  between  the  ApEn  metric  and  the  information  compatibility  (dynamic 
attending)  of  the  presentation  of  data  to  the  operator.  The  ApEn  metric  can  be  calculated  in  real  time 
when  the  dynamic  aspects  of  the  display  may  change.  Both  the  operator  and  experimenter  are  blind  to 
the  active  display  characteristics  and  how  it  evolves  with  time. 


3.  HYPOTHESIS 


This  study  will  show  the  effectiveness  of  the  approximate  entropy  metric  to  identify  when  the  display 
dynamics  may  change.  The  goal  is  to  reject  the  null  hypothesis  Ho  at  an  a  level  of  0.05  such  that: 

Ho:  There  is  no  change  in  the  approximate  entropy  measure  during  extreme  high  and  low  periods  of 
infonnation  incompatibility  (dynamic  attending). 

4.  METHODS 


Data  from  eight  USAF  subjects/contractors  were  collected  as  they  performed  a  time  estimation  task. 
Two  levels  of  information  compatibility  (dynamic  attending)  were  presented  to  the  subjects.  These 
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difficulty  levels  were  kept  blind  from  both  the  experimenter  and  the  subjects.  Figure  1  illustrates  the 
visual  scene  in  which  the  subjects  make  estimates  of  the  time  perception. 


Figure  1  -  Time  Estimation  Task 


In  Figure  1,  the  subject  sees  a  pattern  of  stimuli  being  presented  (small  circle)  in  a  certain  timing 
pattern.  The  small  circle  traverses  the  larger  circular  pattern  (the  large  circle  is  actually  hidden)  until 
the  small  circle  changes  color.  After  this  color  change,  the  subject  has  to  estimate  the  time  for  the  next 
appearance  of  the  circle  as  being  either  “early”,  “late”  or  “on  time”  in  comparison  to  the  prior  timing 
pattern.  The  subject  responds  with  a  mouse  cursor  as  to  one  of  the  three  choices  (early,  late  or  on 
time).  From  the  works  of  Skelly  et  ah,  Table  1  displays  the  four  possible  paradigms  that  are  known  to 
produce  different  levels  of  dynamic  attention  or  discordancy.  Going  down  the  table  is  in  a  direction 
of  increasing  task  difficulty  as  determined  from  prior  studies. 


Figure  2  portrays  an  hypothesized  trial  from  one  possible  set  of  runs.  The  task  may  be  easy  for  a 
certain  period  of  time  and  then,  without  notice  to  the  subject  or  experimenter,  the  tasks  may  enter  a 
region  of  increased  difficulty  for  some  period  of  time.  Finally  the  task  set  may  decrease  its  level  of 
difficulty.  The  opposite  sequence  may  also  occur.  The  presentations  of  the  task  difficulty  levels  are 
counter  balanced  in  the  data  presented  here  to  eliminate  confounding  due  to  fatigue  or  ordering 
effects.  Also  in  Figure  2  is  a  presumed  plot  of  the  variable,  approximate  entropy,  drawn  on  the  same 
time  axis.  The  conjecture  is  that  as  the 


Table  1  -  Four  Types  of  Dynamic  Attending 


Task  Condition 

Difficulty  Level 

Equal  Time 
Equal  Space 

Easiest  Task- 
Level  1 

Equal  Time 

Next  Most  Difficult 
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Unequal  Space 

Task  -  Level  2 

Unequal  Time- 
Unequal  Space 

Third  Most 
Difficult  Task  - 
Level  3 

Unequal  Time  - 
Equal  Space 

Most  Difficult 
Task  -  Level  4 

difficulty  of  the  task  increases,  the  entropy  would  increase,  accordingly.  The  entropy  is  evaluated  by 
the  degree  of  error  made  by  the  subjects  as  they  respond  to  the  question  of  whether  the  last 
appearance  of  the  circle  was  either  early,  late,  or  on  time.  There  are  five  levels  of  error  (zero  error,  + 
one  unit  and  +_two  units  of  error)  that  a  subject  can  make.  The  approximate  entropy  depends  on  the 
absolute  value  of  the  level  of  the  error  made  by  the  subject. 


Task  Approximate 


Figure  2  -  An  Example  of  The  Temporal  Presentation  of  Tasks 

5.  APPARATUS 

A  PC-based  Borlan  C program  operated  in  DOS  was  employed  to  generate  the  experimental  scenario 
and  collect  data.  Data  analysis  was  accomplished  with  MATLAB  ™with  statistical  analysis  conducted 
with  JMP  4.0  (SAS  Institute).  Figure  3  shows  the  facilities  at  the  Air  Force  Research  Laboratory  to  study 
dynamic  attending  issues.  Figure  4  portrays  one  of  the  subjects  in  an  early  pretraining  run. 
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Figure  3  -  Experimental  Setups  to  Perform  Time  Estimation  Testing 


6.  EXPERIMENTAL  DESIGN 

Both  genders  of  healthy  USAF  subjects/contractors  participated  in  this  experiment.  This  research  had  the 
goal  of  correlating  the  approximate  entropy  measure  with  the  task  difficulty  of  the  underlying  dynamical 
experimental  conditions  the  operators  were  being  exposed  to.  One  data  run  of  the  training,  the  easy  or 
hard  task  consisted  of  64  trials.  Frequent  breaks 


Figure  4  -  Subject  Running  During  the  Time  Estimation  Task 


were  given  to  the  subjects  so  that  they  spent  no  more  than  60  minutes  performing  the  testing  each  day. 
The  dependent  measure  is  the  error  (between  the  real  world  and  the  perceived  value)  in  correctly 
classifying  the  differences  in  the  time  estimate  when  the  last  circle  should  appear.  The  approximate 
entropy  measure  was  derived  from  this  error  signal  in  real  time  as  described  in  the  appendix.  Errors 
are  assigned  as  one  unit  if  the  subject  guessed  “On  Time”  when  the  actual  stimulus  was  either  early  or 
late.  An  error  of  two  units  could  be  recorded  if  the  subject  responded  “Late”  when  the  actual  stimulus 
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was  early.  The  reverse  situation  would  also  produce  two  units  of  error,  if  it  occurred.  Other  cases  of  1 
error  unit  were  also  possible. 


7.  RESULTS 


We  report  data  from  eight  subjects  who  have  completed  at  least  one  training  day  and  two  additional 
data  days  of  the  time  estimation  testing.  The  analysis  of  the  data  is  conducted  as  a  within-subjects, 
full  factorial  design.  There  are  two  analyses  to  describe.  The  first  analysis  deals  with  the  empirical 
performance  results  based  on  the  errors  accumulated.  The  second  analysis  deals  with  the  use  of  the 
approximate  entropy  metric  to  distinguish  the  subject’s  response  during  the  presentation  of  three 
levels  of  task  difficulty  (dynamic  attending). 


7. 1  PERFORMANCE  RESULTS 


Similar  to  the  work  reported  earlier  by  Skelly  and  colleagues  (Skelly  (2003),  Skelly  et  al.  (2000)  and 
Jones  and  Skelly  (1995)),  we  describe  some  basic  tests  on  the  data  collected  in  this  study.  Figure  5 
shows  performance  data  from  one  of  the  subjects  during  three  independent  testing  regimes.  The 
absolute  value  of  the  error  signal  is  plotted  on  the  y  axis,  the  trial  number  is  on  the  x  axis. 


Trial  number 

In  Figure  5,  for  the  first  64  trials  (training  run),  the  subject  made  one  error  of  time  judgment  during 
this  event.  This  was  the  final  training  run  and  subjects  were  required  to  have  over  95%  correct  time 
judgments  prior  to  entering  into  the  latter  phases  of  testing.  During  the  training  and  data  runs, 
subjects  had  5  brief  breaks.  The  total  duration  of  a  run  was  19  minutes  to  complete  all  64  trials.  After 
completion  of  the  training  run,  the  subjects  would  either  perform  an  “easy”  task  for  64  trials  or  a 
“hard”  task.  The  distinction  between  easy  or  hard  task  was  defined  by  Dr.  Skelly  (different 
techniques  of  dynamic  attending)  and  was  randomly  presented  to  the  subjects  in  blocks  of  64  trials. 
In  Figure  5,  the  easier  task  was  presented  prior  to  the  more  difficult  task.  In  Figure  5,  the  number  of 
errors  was  found  to  be: 
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Training  Task:  1  error  in  64  trials. 

Easy  Task:  9  errors  in  64  trials. 

Hard  Task:  18  errors  in  64  trials. 

Of  course,  if  a  subject  received  the  treatments  easy  and  then  hard  on  one  data  day,  then  he  must  also 
receive  the  reverse  order  of  hard  and  easy  on  the  next  data  day  to  counterbalance  the  experimental 
conditions.  It  should  be  emphasized  that  the  task  complexity  (bits)  and  the  bit  rate  of  presentation  of 
infonnation  (Fitts’  law  or  bits/second)  were  identical  for  the  easy  and  hard  tasks.  It  was  only  the 
dynamic  presentation  of  the  infonnation  that  was  different  between  the  tasks.  Since  brevity  must  be 
the  style  here,  the  first  analysis  deals  with  performance  in  terms  of  the  number  of  enors  that  occurred. 
This  is  detennined  versus  the  levels  of  task  difficulty  exhibited  in  Table  2.  The  mean  and  SE  of  the 
number  of  errors  across  all  8  subjects  accumulated  during  the  64  trials  is  displayed.  A  one-way 
ANOVA  on  the  difference  of  means  of  the  dependent  measure  of  errors  in  time  perception  across  the 
task  difficulty  levels  is  provided  to  examine  the  efficacy  of  the  experimental  scenario  to  elicit  a 
performance  change. 

Table  2  -  Mean  and  SE  of  Errors  Observed 
Training  Easy  Task  Hard  Task 

Mean  SE  Mean  SE  Mean  SE 

1.91  1.35  9.64  1.35  20.54  1.35 

Table  3  illustrates  the  results  of  a  one-way  ANOVA  (t-test  of  means  -  all  pairs  Tukey-Kramer)  based 
on  JMP  4.0  analysis  (SAS  Institute,  2004). 

Table  3  -  One  Way  ANOVA  for  Errors 
Source  DoF  X  Squares  Mean  Squares  F-Ratio  Prob>F 
Task  Level  2  1928.8  964.394  48.045  <.0001 

Error  30  602.2  20.073 

C.  Total  32  2530.96 

Figure  6  displays  the  output  of  the  comparison  of  means  of  errors  across  all  8  subjects  which  clearly 
(with  the  one-way  ANOVA  results  of  Table  2)  demonstrates  the  efficacy  of  the  experimental 
paradigm  to  elicit  a  performance  change  as  predicted  by  the  prior  works  of  Skelly  et  al. 
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Figure  6  -  Distribution  of  Errors  for  the  Three  Levels  of  Task  Difficulty  as  Defined  by  Skelly  et  al. 


7.2  RESULTS  FROM  ApEn  ANALYSIS 


A  second  goal  of  this  investigation  was  to  show  the  efficacy  of  the  ApEn  metric  to  demonstrate  that 
the  disorder  in  the  operator’s  response  is  also  affected  in  a  dynamic  time  sense.  The  platfonn  for  the 
experimental  design  (dynamic  attending)  discussed  in  section  7.1  clearly  provides  an  investigational 
scenario  to  manipulate  the  disorder  in  the  response  of  the  operator  in  terms  of  task  difficulty.  A  plot 
of  the  ApEn  versus  time  is  shown  in  Figure  7  for  one  subject  and  the  three  runs:  training,  the  easy 
task,  and  the  hard  task.  The  Appendix  describes  the  implementation  issues  for  the  generation  of  the 
ApEn  calculation  portrayed  in  Figure  7.  Other  technical  details  on  the  ApEn  for  the  use  of 
determination  of  irregularity  in  real  time  empirical  data  can  also  be  found  in  Repperger,  et  al.,  2004, 
published  in  this  same  conference  proceedings.  In  Figure  7  it  is  noted  that  three  different  levels  of 
ApEn  seem  to  appear  at  each  level  of  task  difficulty.  Table  4  shows  the  results  of  averages  of  the 
ApEn  metric  over  the  8  subjects.  The  averages  across  subjects  were  conducted  for  the  mean  ApEn 
value  during  the  64  trials  within  the  same  task  difficulty  interval  as  indicated  in  Figure  7. 
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Figure  7  -  Side  view  of  ApEn  Plot  versus  trial  number  for  one  subject. 


Table  4  -  Mean  and  SE  of  ApEn  Calculated 
Training  Easy  Task  Hard  Task 

Mean  SE  Mean  SE  Mean  SE 

0.074  0.02  0.182  0.02  0.247  0.02 


Figure  8  illustrates  a  three  dimensional  plot  of  the  ApEn  function  versus  time  for  the  same  data  as 
displayed  in  Figure  7  with  the  third  axis  representing  r  =  the  threshold  window  as  discussed  in  the 
appendix. 


Figure  8  -  ApEn  Plot  versus  Trial  Number  and  r  =  Threshold  Window 


Table  5  illustrates  the  results  of  the  one-way  ANOVA  (t-test  between  means  -  Tukey-Kramer)  based 
on  JMP  4.0  analysis  (SAS  Institute,  2002)  for  the  ApEn  metric  in  Figures  7-8. 
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Table  5  -  One-Way  ANOVA  for  Magnitude  of  ApEn 


Source  DoF  Z  Squares  Mean  Squares  F-Ratio 


Task  Level  2  0.3031 

Error  30 
C.  Total 


0.1515  21.16 

0.2148  0.0072 

32  0.5179 


Prob>F 

<.0001 


Figure  9  portrays  the  JMP  4.0  analysis  output  as  discussed  in  Table  5  for  a  one-way  ANOVA  on  the 
dependent  measure  of  magnitude  of  ApEn  during  each  of  the  three  levels  of  task  difficulty  considered 
in  this  experimental  study. 


8.  DISCUSSION 


Tables  3  and  5,  presented  in  section  7,  clearly  show  that  the  efficacy  of  the  experimental  paradigm  to 
elicit  a  performance  change  as  well  as  the  ability  of  the  ApEn  metric  to  provide  sufficient  sensitivity 
to  sense  this  change.  Clearly  all  the  levels  of  different  dynamic  attending 


Figure  9  -  Distribution  of  Means  of  ApEn  values  over  levels  of  Task  Difficulty  Across  8  Subjects. 


difficulty  can  be  distinguished  either  through  the  performance  changes  observed  or  via  the  ApEn 
metric  which  can  be  developed  online  in  real  time. 

In  summary,  a  preliminary  study  has  been  conducted  in  time  estimation  when  the  dynamic  display  has 
a  known  change  in  dynamic  attending  capability,  but  with  identical  task  complexity  and  complexity 
rate.  The  efficacy  of  the  approximate  entropy  measure  to  capture  and  characterize  this  dynamic 
change  has  been  evaluated.  This  important  real  time  measure  (ApEn)  provides  a  valuable  means  of 
describing  how  humans  deal  with  complex  systems  which  may  have  characteristics  that  vary  with 
time. 
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9.  CONCLUSIONS 


Studying  dynamic  displays  is  a  difficult  task.  This  preliminary  study  looked  at  a  measure  of  this 
interaction  in  terms  of  the  disorder  of  the  response  of  the  subjects  (entropy  measure).  Additional 
studies  will  examine  display  conditions  and  the  ability  of  this  measure  or  other  metrics  to  capture  the 
true  dynamic  interaction  of  the  operator  within  a  changing  dynamic  environment. 
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APPENDIX  of  [4]  -  DETAILS  REGARDING  IMPLEMENT AION  OF  THE  ApEn  METRIC 

Following  the  other  reference  (Repperger,  et  al.,  2004,  published  in  [3]),  the  choice  was  made  of  m=l, 
n=10,  and  the  window  r  was  adjusted  to  be  approximately  25%  of  the  standard  deviation  of  the  data, 
as  suggested  by  Pincus  and  Viscarello  (1992).  The  difficulty  in  implementing  the  ApEn  measure  is 
because  this  is  a  discrete  task,  rather  than  a  continuous  time  series  sampled  at  a  uniform  sampling 
rate.  To  convert  this  time  estimation  experimental  paradigm  into  a  procedure  amenable  to  analysis 
using  the  ApEn  metric,  one  needs  to  look  at  the  error  signal,  e.g.  as  it  appears  in  Figure  5.  The  kind 
suggestion"  was  made  that  to  convert  a  sequence  of  discrete  tasks  into  a  means  amenable  for 
calculation  by  the  ApEn  metric,  it  is  only  necessary  to  concatenate  the  error  trials  versus  time,  as 
shown  in  Figure  5.  The  ApEn  analysis  considers  the  sequence  of  64  x  3  =  192  trials  as  a  single  time 
series.  The  independent  variable  is  the  trial  number  (1-192).  The  dependent  measure  is  the  absolute 
value  of  the  error  signal  (which  is  0,  1,  or  2  units).  The  ApEn  analysis  was  applied  to  the  data  by 
treating  the  sequence  of  192  trials  as  192  uniform  samples  of  time  data.  Since  the  ApEn  analysis  is 
independent  of  the  sampling  time,  it  does  not  impact  the  analysis  procedure  used  here. 
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Abstract — An  important  problem  in  computer  vision  is  to 
determine  the  orientation  of  a  rigid  body  in  an  image.  This 
can  be  accomplished  by  matching  points  or  line  segments  that 
naturally  appear  on  the  object.  Several  elegant  and  computation¬ 
ally  fast  algorithms  based  on  the  singular  value  decomposition 
and  quaternions  have  been  introduced  to  solve  this  problem. 
In  this  article,  the  authors  first  examine  the  important  special 
case  of  identifying  the  attitude  of  2D  objects  and  introduce  a 
particularly  elegant  solution  based  on  the  mathematical  structure 
of  the  complex  plane.  Motivated  by  this  simple  solution  to  the 
2D  case,  a  new  derivation  of  the  3D  case  based  on  the  polar 
decomposition  is  presented.  This  derivation  is  in  many  ways  more 
natural  than  previous  derivations,  particularly  when  the  model 
and  data  contain  no  noise. 

Index  Terms — Absolute  orientation,  least  squares,  polar  de¬ 
composition. 

I.  INTRODUCTION 

A  fundamental  problem  in  computer  vision  is  the  determina¬ 
tion  of  the  orientation  of  a  rigid  object.  An  effective  approach 
to  this  problem  is  to  match  a  set  of  points  on  the  object  with  the 
corresponding  points  on  a  model.  In  particular,  the  following 
mathematical  problem  appears  in  a  number  of  references  [1]- 
[6]. 

Two  point  sets  {a,}  and  {b, }  of  N  vectors  in  the  plane  or 
in  3-space  are  related  by 

b,  =  Aa,  +  t,  +  n,  (1) 

where  R  is  a  rotation  matrix,  t  is  a  translation  vector,  and  n, 
is  a  noise  term.  The  set  {a;}  corresponds  to  the  location  of 
several  specified  points  on  a  model  of  the  object  while  {b;} 
represents  the  corresponding  points  on  the  object  in  an  image. 
The  goal  is  to  determine  R  and  t  to  minimize 

N 

F(t,R)  =  ^||bt-(f?ai+t)||2  (2) 

*= l 

where  ||  •  ||  is  the  standard  2-norm.  The  vector  t  and  rotation 
matrix  R  represent  the  location  and  orientation  of  the  object  in 
the  image.  The  problem  of  matching  line  segments,  although 
more  complicated,  results  in  essentially  the  same  type  of 
optimization  problem  [2]. 

Several  approaches  to  this  problem  have  been  described 
in  the  literature  including  matrix-based  solutions  such  as  the 


singular  value  decomposition  (SVD)  [l]-[3]  and  quaternion- 
based  solutions  [4], [5],  In  the  next  section,  we  introduce  a 
new  proof  for  the  2D  case  based  on  simple  properties  of 
complex  numbers.  In  particular,  it  is  shown  that  the  solution 
for  the  orientation  of  the  object  is  given  by  the  polar  form 
of  a  particular  complex  number.  In  Section  III,  this  solution 
motivates  a  natural  solution  to  the  general  3D  problem  based 
on  the  polar  decomposition  of  a  particular  matrix.  In  fact,  the 
polar  decomposition  is  the  obvious  solution  when  no  noise  is 
present  in  the  problem.  After  solving  for  the  case  when  there 
is  noise  in  the  data,  we  provide  a  new  proof  that  the  same 
solution  holds  when  there  is  noise  in  both  the  model  and  the 
data.  Lastly,  conclusions  appear  in  Section  IV. 

II.  SOLVING  THE  PLANAR  CASE  USING  THE 
COMPLEX  PLANE 

A  popular  solution  to  the  orientation  problem  is  based  on 
the  singular  value  decomposition  (SVD).  However,  as  one 
would  expect,  the  much  simpler  2D  case  does  not  require 
the  sophistication  of  an  SVD,  not  only  because  of  the  smaller 
dimension  size,  but  more  importantly,  because  of  the  commu¬ 
tativity  of  the  rotation  operation.  In  this  case,  it  is  convenient  to 
formulate  the  problem  in  terms  of  complex  numbers.  Suppose 
that  a  =  [  ax  ay  ]  and  b  =  [  bx  bv  ]  are  vectors 
in  the  plane.  If  we  write  these  vectors  in  complex  number 
notation  as  a  =  ax  +  jay  and  b  =  bx  +  jby,  then  the  inner 
product  a  •  b  of  the  two  vectors  is  given  in  complex  number 
notation  as  R.e(a*b)  where  a*  denotes  the  complex  conjugate 
of  the  complex  number  a  and  where  Re(z)  denotes  the  real 
part  of  2.  Furthermore,  the  norm  squared  ||a||2  of  the  vector  a 
is  given  by  |a|2  =  a* a,  and  the  rotation  matrix  corresponding 
to  a  counterclockwise  rotation  of  9  radians  is  given  by  the 
complex  number  e-7  .  Based  on  this  formulation,  the  goal  is  to 
minimize  the  objective  function 

1  N 

F(t,e^e)=-Y/\bi-(ejeai  +  t)\2  (3) 

*= l 

where  the  complex  number  t  and  the  real  number  0  S  [0.  27t) 
represent  the  position  and  orientation  of  the  object,  respec¬ 
tively.  Since  determining  the  optimal  to(9)  for  a  given  9  is 
a  routine  least  squares  calculation,  we  merely  state  the  result 
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that 

t0  (6>)  =b-ej8a  (4) 

where  b  =  jj  J^iLi  k  and  a  =  J2tLi a* •  We  define 

Hi  =  a-i  —  a  and  bi  =  bi  —  b  and  say  that  {a.j}  and  {6;} 
are  the  unbiased  versions  of  {a,}  and  {6,},  respectively.  It 
then  follows  that 


F{ffe)  =  fc-e^l2 


N 

1  sr^N 

~  N  Z^i=l 


Since  ±  Y^=i 


maximize  Re 


\bi\ 2  +  |a» 


\h\2  +  |d,:|2  -  2Re((bi)*ejeai) 

J(5) 

is  fixed  by  the  data,  we  want  to 


jvEi=i(M*«i  e' 


is  clearly  achieved  by  choosing 


J8 


with  respect  to  6,  which 


N 


0  =  -  arg  ^2(bt)*at  )  =  arg  ^  5*b 


N 


(6) 


where  arg(z)  denotes  the  argument  of  the  complex  number 
z  and  has  range  [0,  27t).  The  location  and  orientation  of  the 
object  is  then  given  by  (4)  and  (6),  respectively. 

Another  approach,  that  will  serve  as  a  guide  in  the  next 
section  to  solve  the  3D  case,  is  to  write 


F{t,e?e) 


1 

N  2-^i= 1 

N2  +  N2 

1  y^AT 

N  2^i= 1 

N2  +  |di|2 

2R  e(c*ej0) 

|c| 2  —  1  +  \c  —  e^l'2 


(7) 

where  c  =  y=1  a* bi.  We  thus  want  to  minimize  | c  —  e J 0 1 2 , 

which  is  clearly  achieved  by  9  =  arg (c)  or,  equivalently,  9  = 

arg  (  E*Li  • 


III.  SOLVING  THE  GENERAL  CASE  USING  THE 
POLAR  DECOMPOSITION 

A  popular  approach  to  solving  the  3D  case  is  based  on 
quaternions  [4],  [5],  While  quaternions  are  a  generalization 
of  complex  numbers,  the  complex  number  approach  of  the 
previous  section  more  naturally  leads  to  a  matrix  solution  based 
on  the  polar  decomposition.  We  begin  by  reformulating  the 
problem  in  matrix  notation  by  letting  A  =  [  ai  •  •  •  a  y  ] , 
B  =  [  bi  •  •  •  bjv  ],  and  N  =  [  ni  •  •  •  njy  ].  Equa¬ 
tion  (1)  then  becomes 

B  =  RA  +  teT  +  N  (8) 

where  e  =  [  1  •••  1  ]T.  The  optimization  problem  then 

becomes  to  minimize 

F(t,  R)  =  \\B  —  (RA  +  teT)|||.  (9) 

subject  to  R  being  a  rotation  matrix  where  ||  •  ||  F  denotes  the 
Frobenius  norm,  which  is  given  by  the  square  root  of  the  sum 
of  the  squares  of  the  matrix  elements. 


A.  The  Noise-free  Case 

We  first  examine  the  simplest  possible  case,  i.e.,  when  no 
noise  is  present.  In  this  ideal  case,  we  have  an  exact  equality, 
which  can  be  written  in  matrix  form  as 

RA  +  teT  =  B  (10) 

where  e  =  [  1  •••  1  ]T.  The  translation  term  t  can  be 

found  by  post-multiplying  (10)  by  7e  to  obtain 

t  =  b  —  R&  (11) 

where  a  =  and  b  =  Writing  A  = 

A— ae2  and  B  =  l>  ber,  we  have  RA  =  B  or,  equivalently, 
AT Rt  =  if  .  Lastly,  pre-multiplying  both  sides  by  A  gives 

AAtRt  =  ABt.  (12) 

We  assume  that  the  3x3  matrix  ABT  has  full  rank.  Then  the 
unique  polar  decomposition  of  ABT  =  PQ  is  given  by  the 
left  hand  side  of  (12)  where  P  =  AAT  is  positive  definite  and 
Q  =  if  is  an  orthogonal  matrix.  Since  P  =  ,4,4 7  is  positive 
definite,  it  follows  that  the  determinants  of  AB1  and  Q  have 
the  same  sign  so  that  R  =  QT  is  a  rotation  matrix  if  and  only 
if  det (Abt)  >  0.  If  det (Abt)  <  0,  then  R  is  a  reflection 
matrix  and  the  orientation  problem  is  ill-defined. 

The  polar  decomposition  is  a  natural  solution  to  the  problem 
when  no  noise  is  present.  Before  continuing  to  the  more 
general  case,  we  observe  that  the  solution  is  particularly  simple 
when  A  has  the  property  that  AAT  =  kl.  In  that  case,  we 
merely  scale  ABT  to  obtain  an  orthogonal  matrix.  An  example 
of  this  occurs  when  the  columns  of  A  correspond  to  the  vertices 
of  a  platonic  solid.  Unfortunately,  choosing  a  data  matrix  so 
that  A  has  this  property  may  result  in  a  numerically  unstable 
solution  when  noise  is  present. 

B.  Noise  in  Only  the  Data 

We  now  return  to  the  problem  formulated  in  (9),  which  can 
be  rewritten  as 

F(t,R)=\\teT  ~(B-RA)\\2f.  (13) 

The  optimal  solution  for  t  for  a  fixed  rotation  matrix  R  is 
given  by  the  pseudoinverse  solution 

t  =  (B  —  RA)  (eT)+  =  ^  (B  -  RA)e  (14) 

so  that  it  follows  that  t  is  once  again  given  by  (11).  Substituting 
this  expression  into  (13)  gives 

\\RA-B\\2f  =  \\ArF-2tr(BTRA)  +  \\BrF 

=  \\A\\2F  +  \\B\\l~2tr(ABTR).  > 

Applying  the  slip-in/slip-out  method,  we  can  write  our  objec¬ 
tive  function  as 

\\ra-b\\2f  =  \\A\\%+m2F-\\ABT\\2F-\\m2F 

+  \\ABT\\2F-2tr(ABTR)+\\R\\2F 
=  \\A\\2F+\\B\\2F-\\ABT\\2F-n 
+\\Abt  —  R\\F. 

(16) 
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Note  the  similarity  of  (16)  with  (7).  The  first  four  terms  in 
the  last  equality  of  (16)  are  independent  of  R,  so  our  problem 
becomes  the  optimization  of  j|H.BT  —  f?||F. 

Note  the  special  structure  of  this  final  version  of  the 
optimization  problem.  We  seek  an  orthogonal  matrix  that  is 
closest  to  the  3x3  matrix  ABT .  More  generally,  we  consider 
the  problem  of  determining 


Uq  =  arg  min  \\M  —  U\\p  (17) 

UeO(n) 

where  M  is  an  a  x  n  full  rank  matrix.  We  solve  this  problem 
by  first  examining  two  special  cases.  First,  suppose  M  is  a 
diagonal  matrix:  M  =  diag(di, . . . ,  dn).  Then 

\\M-U\\%  =  \\D-U\\l 

=  \\DfF  +  \\U\\2F-2tv(DU) 

n  n 

=  y ^(%  +  n  -  2^ ^djUg.  (18) 

i— 1  i— 1 

Since  di  an<3  n  are  fixed,  we  need  to  maximize 

1  diUu.  Since  —1  <  u.a  <  1,  we  have  that 

Z)"=i  diua  <  YJi= i  Mil  and  it  follows  that  Uq  = 

diag(sgn(di), . . .  ,sgn(dn))  G  0(n).  Furthermore,  if  D  is 
positive  definite,  i.e.,  if  di  >  0  for  i  =  1, . . . ,  n,  then  Uq  =  I. 

For  our  second  case,  suppose  that  M  is  symmetric.  Then 
we  can  write  M  as  M  =  VDVT  where  V  G  0(n)  and  D  = 
diag(di, . . . ,  dn )■  We  then  have 

|| M  -  U\\F  =  ||  VDVt  -  U\\F  =  || D  -  VUVt\\f  (19) 

where  the  second  equality  in  (19)  follows  from  the  fact 
that  the  Frobenius  norm  is  invariant  under  pre-  and  post¬ 
multiplication  by  orthogonal  matrices.  By  the  first  case,  it 
follows  that  the  optimal  Uo  G  0(n)  is  given  by  VtUqV  = 
diag(sgn(d1), . . . ,  sgn(d„)),  i.e.,  the  optimal  orthogonal  ma¬ 
trix  is  given  by  Uo  =  Fdiag(sgn(di), . . . ,  sgn(d„))FT.  For 
the  important  case  when  M  is  symmetric  positive  definite,  this 
becomes  Uq  =  I. 

We  are  now  ready  for  the  general  case  when  M  is  an 
arbitrary  full  rank  nx  n  matrix.  In  this  case,  we  can  write  M 
in  its  polar  form:  M  =  PQ  where  P  is  symmetric,  positive 
definite  and  Q  G  0(n).  Then  || M  —  U\\f  =  ||PQ  —  U\\f  = 
|| P  —  UQt\\f*  which  is  minimized  over  U  G  0(n)  by  the 
orthogonal  matrix  Uo  where  UqQt  =  /,  i.e.,  Uq  =  Q- 

We  thus  conclude  that  the  optimal  solution  for  determining 
the  orientation  of  the  object  in  question  is  given  by  the 
orthogonal  matrix  R  =  QT  where  PQ  is  the  unique  polar 
decomposition  of  ABT .  We  note  that  when  the  polar  decom¬ 
position  has  been  mentioned  in  the  literature  as  a  solution, 
it  has  primarily  appeared  as  an  afterthought  of  the  SVD 
solution.  While  these  solutions  are  arguably  equivalent,  the 
SVD  solution  does  not  appear  to  be  as  natural  as  the  polar 
decomposition  solution.  It  is  important  to  note  once  again  that 
R  is  a  rotation  matrix  if  and  only  if  det(Q)  =  1;  otherwise, 
R  is  a  reflection  matrix. 


C.  Noise  in  Both  the  Model  and  Data 


We  next  examine  the  case  when  there  is  noise  in  both  the 
model  and  the  data.  In  particular,  we  consider  the  following 
problem  described  in  [6],  given  here  with  slightly  different 
notation.  Suppose  we  have  two  sets  of  N  noisy  observations 
given  by  two  3  x  N  matrices  A  and  B.  We  assume  that  the 
correspondence  problem  has  already  been  solved  so  that  the 
corresponding  points  are  in  the  same  order  in  A  and  B  and  that 
the  translation  of  the  object  has  already  been  determined.  Our 
problem  then  is  to  find  a  rotation  matrix  R  and  perturbations 
5A  and  6B  which  satisfy 


R{A  +  SA)  =  B  +  SB  (20) 


such  that  ||cL4||!n-|-||<5.B||!>  is  minimized.  The  perturbations  5 A 
and  SB  correspond  to  noise  in  the  model  and  data,  respectively. 
This  problem,  presented  in  slightly  different  notation,  was 
solved  by  Goryn  and  Hein  in  [6],  but  the  solution  presented 
there  relies  heavily  on  the  introduction  of  some  non-obvious 
substitutions.  We  present  a  more  natural  and  intuitive  deriva¬ 
tion.  We  begin  by  first  noting  that  unlike  SB,  the  term  SA  ap¬ 
pears  in  (20)  as  RSA,  suggesting  that  it  may  be  better  to  rewrite 
the  cost  function  as  ||<L4||J,  +  ||5B|||,  =  ||i?<5H|||,  +  ||<5i?|||,, 
where  equality  follows  from  the  fact  that  R  is  orthogonal.  Next, 
equation  (20)  can  be  written  as  RSA  —  SB  =  —  ( RA  —  B ). 
We  thus  want  to  minimize  ||i?i5H|||.  +  ||<5.B||F  subject  to 
RSA-  SB  =  -(RA-  B). 

The  scalar  version  of  the  preceding  problem  statement 
suggests  a  solution  to  this  constrained  optimization  problem. 
In  the  scalar  case,  we  want  to  minimize  x2  +  y 2  over  the 
scalars  x,y,  z  subject  to  the  constraint  ax  +  by  =  f(z)  where 
f(z)  =  — cz  +  d  and  where  a,  b ,  c,  d  are  fixed  parameters. 
Geometrically,  for  a  fixed  z,  the  constraint  can  be  interpreted  as 
the  equation  of  a  line  and  the  optimal  solution  over  x,  y  would 
then  correspond  to  the  point  on  that  line  which  is  closest  to 
the  origin.  We  then  want  to  choose  z  to  place  the  line  as  close 
to  the  origin  as  possible,  which  is  achieved  by  minimizing 
\f(z)\  =  \cz  —  d\.  Once  this  is  done,  the  optimal  x  and  y  can 
be  determined.  This  suggests  minimizing  the  norm  of  RA  —  B 
in  the  non-scalar  case. 

We  now  present  a  formal  derivation  of  the  solution  for  the 
non-scalar  case.  The  constraint  (20)  in  the  general  case  can  be 
written  in  matrix  notation  as 


J _/ 

zAt  ' 

'  RSA  ' 

L  V21 

V21  \ 

SB 

tAra~B) 


(21) 


where  the  l/\/2  term  is  included  so  that  the  rows  of  the 
matrix  on  the  left  are  not  only  mutually  orthogonal,  but  are 
also  normalized.  This  suggests  augmenting  the  matrix  so  that 
it  becomes  orthogonal: 


1 

'  i 

-i ' 

'  RSA  ' 

1 

'  -(RA-B)  ' 

71 

i 

i 

SB 

“71 

RSA  +  SB 

Since  the  Frobenius  norm  is  invariant  under  multiplication  by 
orthogonal  matrices,  we  have  that  the  cost  function  ||(5H|||,  + 
||(5-B||f.  is  given  by 

\\R6A\\2f  +  \\6B\\2f=  ±\\RA-BfF+±\\R5A  +  6B\\2F,  (23) 
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which  is  clearly  minimized  by  setting  SB  =  —RSA  and 
minimizing  ||/L4— B\\2F  over  the  family  of  orthogonal  matrices 
R  so  that  R  is  determined  in  the  same  manner  as  before. 

IV.  CONCLUSION 

In  this  article,  we  have  presented  the  polar  decomposition 
as  the  natural  method  for  solving  the  absolute  orientation 
problem.  Previously,  the  polar  decomposition  approach  was 
mentioned  in  the  literature  as  merely  an  afterthought  of  the 
SVD  approach.  The  polar  decomposition  method  for  the  3D 
case  was  motivated  by  a  complex  number  approach  to  the 
2D  case,  which  is  interesting  in  its  own  right.  Lastly,  we 
have  provided  a  simpler  and  more  natural  proof  that  the  same 
solution  also  holds  for  a  noisy  model. 

Acknowledgment 

The  authors  gratefully  acknowledge  the  support  of  the  Air 
Force  Research  Lab. 


References 

[1]  K.  S.  Arun,  T.  S.  Huang,  and  S.  D.  Blostein,  “Least- squares  fitting  of 
two  3-D  point  sets,”  IEEE  Trans.  PAMI,  vol.  9,  no.  5,  pp.  698-700,  Sept. 

1987. 

[2]  B.  Kamgar- Paris  and  B.  Kamgar-Paris,  “Matching  sets  of  3D  line 
segments  with  application  to  polygonal  arc  matching,”  IEEE  Trans. 
PAMI,  vol.  19,  no.  10,  pp.  1090-1099,  Oct.  1997. 

[3]  B.  K.  P.  Horn,  H.  M.  Hilden,  and  S.  Negahdaripour,  “Closed-form 
solution  of  absolute  orientation  using  orthonormal  matrices,”  Journal 
of  the  Optical  Society  of  America  A,  vol.  5,  no.  7,  pp.  1127-1638,  July 

1988. 

[4]  B.  K.  P.  Horn,  “Closed-form  solution  of  absolute  orientation  using  unit 
quaternions,”  Journal  of  the  Optical  Society  of  America  A,  4(4):629-642, 
April  1987. 

[5]  M.  D.  Wheeler  and  K.  Ikeuchi,  “Iterative  estimation  of  rotation  and 
translation  using  the  quaternion,”  Technical  Report  CMU-CS-95-215, 
Carnegie-Mellon  University,  December  1995. 

[6]  D.  Goryn  and  S.  Hein,  “On  the  estimation  of  rigid  body  rotation  from 
noisy  data,  ”  IEEE  Trans.  PAMI,  vol.  17,  no.  12,  pp.  1219-1220,  Dec. 
1995. 


Approved  for  public  releaS§;  distribution  is  unlimited 


